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Massive gravity is a good theoretical laboratory to study modifications of General Relativity. The 
theory offers a concrete set-up to study models of dark energy, since it admits cosmological self- 
accelerating solutions in the vacuum, in which the size of the acceleration depends on the graviton 
fT^ ' mass. Moreover, non-linear gravitational self-interactions, in the proximity of a matter source, 

manage to mimic the predictions of linearised General Relativity, hence agreeing with solar-system 
("^ ' precision measurements. In this article, we review our work in the subject, classifying, on one 

^SJ ' hand, static solutions, and on the other hand, self-accelerating backgrounds. For what respects 

; I , static solutions we exhibit black hole configurations, together with other solutions that recover 

rS ,■ General Relativity near a source via the Vainshtein mechanism. For the self-accelerating solutions 

we describe a wide class of cosmological backgrounds, including an analysis of their stability. 
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I. INTRODUCTION 

O ' Can the graviton have a mass? A graviton mass breaks the diffeomorphism invariance of General Relativity (GR), 

but has the advantage to potentially provide a theory of dark energy that explains, in a technically natural way, the 
present day acceleration of our Universe. At large scales, gravity is modified with respect to GR, and the theory 
admits cosmological accelerating solutions in which the size of acceleration depends on the graviton mass. This way 
to explain cosmological acceleration is technically natural in the 't Hooft sense, because in the limit of graviton mass 
going to zero one recovers the full diffeomorphism invariance of GR: hence, corrections to the size of dark energy must 
be proportional to the (tiny) graviton mass itself. 

Not aware of all possible consequences of massive gravity for cosmology, Fierz and Pauli (FP), back in 1939, started 

the theoretical study of massive gravity from a field theory perspective [1]. They considered a mass term for linear 

'^ • gravitational perturbations, which is uniquely determined by requiring the absence of ghost degrees of freedom. The 

^T^ [ mass term breaks the gauge (diffeomorphism) invariance of GR, leading to a classical graviton with five degrees of 

freedom, instead of the two found in GR. There have been intensive studies into what happens beyond the linearized 

theory of FP. In 1972, Boulware and Deser (BD) found a scalar ghost mode at the nonlinear level, the so called sixth 

degree of freedom in the FP theory [2] . This issue has been re-examined using an effective field theory approach, where 

gauge invariance is restored by introducing Stiickelberg fields [3]. In this language, the Stiickelberg fields physically 

\^ , play the role of the additional scalar and vector graviton polarizations. They acquire nonlinear interactions which 

Cy_' contain more than two time derivatives, signaling the existence of a ghost [3]. In order to construct a consistent 

theory, nonlinear terms should be added to the FP model, which are tuned to remove the ghost order by order in 

perturbation theory. Interestingly, this approach sheds light on another famous problem with FP massive gravity; due 

to contributions of the scalar degree of freedom, solutions in the FP model do not continuously connect to solutions 

in GR, even in the limit of zero graviton mass. This is known as the van Dam, Veltman, and Zakharov (vDVZ) 

discontinuity [4, 5]. Observations such as light bending in the solar system would exclude the FP theory, no matter 

how small the graviton mass is. In 1972, Vainshtein [6] proposed a mechanism to avoid this conclusion; in the small 

mass limit, the scalar degree of freedom becomes strongly coupled and the linearized FP theory is no longer reliable. 

In this regime, higher order interactions, which are introduced to remove the ghost degree of freedom, should shield 

the scalar interaction and recover GR on sufficiently small scales. 

Until recently, it was thought to be impossible to construct a ghost-free theory for massive gravity that is compatible 
with current observations [7, 8]. Using an effective field theory approach, one can show that in order to avoid the 
presence of a ghost, interactions have to be chosen in such a way that the equations of motion for the scalar and vector 
component of the Stiickelberg field contains no more than two time derivatives. Recently, it was shown that there is 
a finite number of derivative interactions for scalar lagrangians that give rise to second order differential equations. 
These are dubbed Galileon terms because of a symmetry under a constant shift of the scalar field derivative [9]. 
Therefore, one expects that any consistent nonlinear completion of FP contains these Galileon terms, at least in an 
appropriate range of scales in which the scalar dynamics can be somehow isolated from the remaining degrees of 



freedom; this is the so-called decoupling limit [3] . This turns out to be a powerful criterium for building higher order 
interactions with the desired properties. Indeed, following this route, de Rham, Gabadadze and ToUey constructed 
a family of ghost-free extensions to the FP theory, which reduce to the Galileon terms in the decoupling limit. We 
refer to the resulting theory as A3 massive gravity [10]. 

In this article, we review our work to build and analyze exact solutions in A3 massive gravity. As we have briefly 
explained, non-linear effects play an essential role to characterize phenomenological consequences of this theory. Then, 
the analysis of exact solutions of the equations of motion, obtained by imposing appropriate symmetries (spherical 
symmetry for static space-times, or homogeneity and isotropy for cosmological set-ups), make manifest, in idealized 
but representative situations, how the non-linear dynamics of the graviton degrees of freedom respond to the presence 
of a source, or, at very large scales, to the graviton mass itself. After all, looking back to the past, we know that 
the knowledge of exact solutions of non-linear field equations have been of crucial importance to understand GR. 
The Schwarzschild solution lead to the discovery of the concept of black hole, and play an essential role for analyzing 
the dynamics of objects around massive sources in GR; and modern cosmology would be unthinkable without the 
use of Friedmann-Robertson- Walker solutions. Exact solutions in massive gravity might lead to the discovery and 
understanding of new features and concepts in a theory of gravitation that can lead to important developments for 
our comprehension of gravitational interactions. 

This review is organized as follows: in Section II we explicitly construct the A3 massive gravity theory, while in Sec- 
tion III the most general Ansatz for spherically symmetric solution is introduced. This leads to two branches of static 
solutions: one exhibiting the Vainshtein mechanism, and the other representing a generalisation of Schwarzschild- 
(A)dS black holes. In Section IV we explore cosmological self-accelerating solutions and their stability under pertur- 
bations. Finally, in Section V we conclude, also outlining possible directions for future research. 

II. GHOST-FREE MASSIVE GRAVITY 

We begin with the covariant Fierz-Pauli mass term in four-dimensions, given by 

Cpp^m'^—gU^^K U^^^^{H,,H^'^-H^), (1) 

where to is a parameter with units of mass and the tensor iJ^^ is a covariantisation of the metric perturbations, 
namely 

5^1/ = ??A"^ + ^M^ = ^A"^ + S^t., with Y.f^i, = df^(j)"d^cj)'\ai3- (2) 

The Stiickelberg fields 0" are introduced to restore reparametrisation invariance, hence transforming as scalar from 
the point of view of the physical metric [3]. The internal metric rjap corresponds to a non-dynamical reference metric, 
usually assumed to be Minkowski space-time. The dynamics of the Stuckelberg fields (/>" are at the origin of the two 
features discussed in the introduction: the BD ghost excitation and the vDVZ discontinuity. With respect to the first 
issue, as noticed by Fierz and Pauli, one can remove the ghost excitation, to linear order in perturbations, by choosing 
the quadratic structure iJ^^iJ^"^ — H^ . When expressed in the Stiickelberg field language, terms in the action are 
arranged in a such way to constraint one of the four Stuckelberg fields to be non-dynamical. However, when going 
beyond linear order, this constraint disappears, signaling the emergence of an additional ghost mode [3]. Remarkably, 
Ref. [11] has shown how to construct a potential, tuned at each order in powers of H^^, to hold the constraint and 
remove one of the Stiickelberg fields. Even though the potential is expressed in terms of an infinite series of terms for 
H^j_i,, it can be resummed into the following finite form [10, 12] 

U = — TO,^ [W2 + 0^3 W3 + 0^4 U4\ , (3) 

where q;„ are free dimensionless parameters, Lin — n\ det„(/C) and the tensor /C "^ is defined as 

/c/ ^ c - ( VF^) " . (4) 



(The square root is formally understood as v A^u v ^a — ^u ■) The relationship of these potentials with a determinant 
resides on the following property, which holds for squared real matrices and a complex number z 

00 
det(I+z/C) = l + y z"det(/C), (5) 



where each determinant can be written in terms of traces as 
det(/C) = tr/C, 
det(/C) = (tr/C)2 - tr(/C2), 
det(/C) = (tr/C)3-3(tr/C)(tr/C2)+2tr/C^ 
det(/C) = (tr/C)^ - 6(tr/C)2(tr/C2) + 8(tr/C)(tr/C^) + 3(tr/C2)2 - Gtr/C* . 

All terms det„(/C) with n > 4 vanish in four dimensions. If one chooses a sum of determinants of the form J2n=i det(I+ 
ZnIC) — 4, one can generate each det„(/C) term with a separate coefficient a„, provided a solution to J2i=i ^n — '^n 
exists, which is guaranteed by the Newton identities. Therefore, the massive gravity theory can be written in full as 

C=^^./—g{R-2A-m^U), (6) 

where U is given by (3) and we have introduced an additional bare cosmological constant A. Notice that in order to 
obtain the Fierz Pauli term (1) as the first order correction to GR, we have ignored the tadpole term deti(/C). 

We can write the equations of motion in a more familiar way by using the potential term (3) as the source of 
peculiar energy momentum tensor. In this way, the Einstein equations read 



^^v — -f^jy; \j) 



where the energy momentum tensor is defined as 






The theory defined by (6) has Minkowski spacetime as trivial solution when A = 0, hence one can rewrite the metric 
g^i, and the scalars (^^ as deviations from flat space, namely 

gt^v^ViJ^v + h^.u, 4>" ^ x" - tt"- , (9) 

where x" are the usual cartesian coordinates spanning rjafj. In what follows we will use cff^ or tt", having in mind 
that (9) relates them. Moreover, the unitary gauge (where tt^ = 0) simplifies the potential (3) considerably, and we 
will start with this choice in what follows. 

There have been intensive studies in the issue of BD ghost in this theory [13]. The general (but not universal [14]) 
consensus is that there is indeed no BD ghost and the maximum number of propagation modes in this theory is five. 
However, this does not preclude a possibility that one of the five modes becomes a ghost around some backgrounds. 

III. SPHERICALLY SYMMETRIC SOLUTIONS 

In this section, we review spherically symmetric solutions in the unitary gauge by following Refs. [12, 15]. The most 
general Ansatz with spherical symmetry, before fixing the gauge, is 

ds^ = -b{t, r)2 df + a{t, rf dr^ + 2d(t, r) dtdr + c(i, rfdV?, (10) 

where d'^ = d9^ + sin OdcjP , and the Stuckelberg fields have the structure 

</)0 = /(t,r), ^^^g(t,r)-. (11) 

r 

We start focussing on the unitary gauge (tt'' = or, equivalently, f ~ t and g ~ r) and look for static solutions that 
do not depend, explicitly, on time. The metric ansatz (10) reduces to 

ds^ = -b{rf df + a{rf dr^ + 2d{r) dtdr + cirfdO^ . (12) 

Furthermore, we choose to write the non-dynamical flat metric in (2) as ds^ = —dt^ + dr^ + r'^d^^ . It should be 
noticed that this is not a coordinate choice, but a way to simplify the expressions. Indeed, we have chosen the unitary 
gauge, hence we are left with a theory which is not diffeomorphism invariant. Hence, in this context physics does 



depend on the choice of coordinates: we wih further explore this fact in section IV. Any change of coordinate normaUy 
breaks the unitary gauge and switches on a non-trivial profile for the Stiickelberg fields. This also implies that the 
static solutions considered in the unitary gauge do not provide all the spherically symmetric and static solutions in 
this theory. Other static solutions might exist with non-trivial Stiickelberg fields turned on. 

Conscious of these limitations, let us start with this gauge choice: we will relax it in what follows. We plug 
the previous metric into the Einstein equations (7), and observe that the Einstein tensor G^i, satisfies the identity 
d{r) Gtt + b{r)'^ Gtr — 0, which implies the algebraic constraint — d{r) Tjj + 6(r)^ Tj^. This last equation implies 

rf(r)(cor-c(r)) - 0, (13) 

where Cg is a function of a^ and 0:4 only (see Section IIIB). This constraint is solved in two possible ways, defining 
two branches of solutions: 

d{r) ~ or c{r) = c^r . (14) 

In the next sections we will analyze each of these two branches separately. We will start from the diagonal 
one d{r) = in Section III A, where the Vainshtein eff'ect takes place and can be analyzed in a systematic way. 
Then we will proceed in Section III B to study the class of solutions with a non-diagonal metric and c{r) = cq r, 
corresponding to non-asymptotically flat, Schwarzschild-(Anti)-de Sitter solutions that can be relevant to explain 
present-day cosmological acceleration. 

A. Branch I: Vainshtein mechanism at work 

The problem of finding exact vacuum solutions in this branch d{r) = is an open question, but we can make 
interesting progresses by considering perturbations (not necessarily small) from flat space. The following Ansatz is 
useful 

b{r) = l + N{r)dt^, a{r) = {l + F{r)y^''^ c{r) = {1 + H{r)y\ (15) 

Furthermore, it is convenient to introduce a new radial coordinate p = Y+57rT' ^'~' ^^^^ the linearised metric is 
expressed as 

ds^ = -(1 + n)dt^ + (1 - /) V + p^dfl^, (16) 

where f{p) — F(r(/9)) — 2h(p) — 2ph'{p), n{p) — 2iV(r(p)), h{p) = H(r{p)^ and a prime denotes a derivative with 
respect to p. As discussed above, one should be careful with this change of coordinates since, after fixing a gauge, a 
change of frame in the metric breaks the unitary gauge and switches on the Stiickelberg fields tt'' . It turns out that 
this coordinate transformation excites the radial component of tt^, which explicitly reads tt^ = ph. Therefore, from 
now on one can think of h as geometrically corresponding to the only non-zero component of the Stiickelberg field 
TT^. At linear order, the equations for the functions n{p), f{p) and h{p) in the new variable radial p are 

= {m^p'^ + 2) f + 2p{f' + m^p^h' + 3m^ph) , (17) 

= ^m^p^{n-4:h)- pn' - f, (18) 

= f + ^pn'. (19) 

In this linear expansion, the solutions for n and / are 

ACM 
f = -^— (l+mp)e-"^ (21) 

3p 

where we fix the integration constant so that M is the mass of a point particle at the origin, and 8itG = MZ . These 
solutions exhibit the vDVZ discontinuity, since the post-Newtonian parameter 7 = f/n is 7 = ^(1 -|- mp), which in 
the massless limit reduces to 7 = 1/2, in disagreement with GR predictions and solar system observations. 

However, in order to understand what really happens in this limit, we must carefully analyse the behaviour of the 
function h(r) as m — >■ 0. The study of the equations of motion for this metric component makes manifest the non- 
linear effects responsible for the Vainshtein mechanism. To do this, we consider scales below the Compton wavelength 



mp <C 1, and at the same time ignore higher order terms in GM . Under these approximations, the equations of motion 
can still be truncated to linear order in / and n, but since h is not necessarily small, we have to keep all non-linear 
terms in h. In other words, we take the usual weak field limit for the metric fields, but keep all non-linearities in the 
component h, since we expect regions where non-linear effects in h become important. As shown in [12], the field 
equations reduce to the following system of coupled equations for the fields /, n, h: 
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Equation (24) is a quintic algebraic equation in h, except for the special case where (3 — 0, where it reduces to a 
cubic equation. Thus, after obtaining a solution for h from equation (24), one can calculate the gravitational potentials 
/ and n using (22) and (23). In the particular case of /3 = 0, it is possible to describe the solution in a simple way 
[15]. For large radial p values, one can linearise the equations in h, recovering the solution in Eqs. (20)-(21), to 
first order in mp. On the other hand, the Vainshtein mechanism applies, and below the so-called Vainshtein radius, 
Pv = {GMm~'^Y'^, h becomes larger than one and the non-linear terms in h in cq. (24) become important, recovering 
GR close to a matter source. Actually, for p -^ pv the solution for h is simply given by \h\ — pv /{ap) ^ 1. The 
latter solution for h and Eq. (24) with /3 = implies 



2pn' 
I 



2GM 



1 + - A ■ 
2a \pvj 




(29) 



Therefore, corrections to the GR solutions are indeed small for p < py, as shown in the left plot of Fig. 1. Note 
that if we consider a finite size matter source, it was shown that there is no stable solution that interpolates from the 
Vainshtein region to the asymptotically flat solution, and the Vainshtein region is naturally matched onto a solution 
that asymptotes to a non-flat cosmological background [16]. 
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FIG. 1: Numerical solution for drf = /', drti — n' , and the quotient 7' = f' /2n' around the Vainshtein radius p„ (left) and 
the Compton Wavelength p ~ 1/m (right). Region 1 (2) shows how GR solutions are (not) recovered inside (outside) the 
Vainshtein radius pv- Region 3 shows the asymptotic decay of the linear solutions (Eqs. (20)-(21)). Here, GM = 1 and a — 1. 



For the /3 7^ case, analytic solutions of the previous algebraic equations can not be found in general. In the 
following we will follow the approach presented by Ref. [17]. It is possible to determine exactly how many local 
solutions exist in a neighbourhood of infinity at p = +00, which we refer as asymptotic solutions, and moreover how 



many local solutions exist in a neighbourhood of p ~ 0^ , which we call inner solutions. Furthermore, we can find 
analytically their leading behaviour as a function of p. Any global solution of (24) should necessarily interpolate 
between one of the asymptotic solutions and one of the inner solutions. Therefore, our aim is to understand, for each 
point in the (a, j3) phase space, whether and how the above solutions match. 

In a neighbourhood of p — > +cxi there are, depending on the value of (a, /3), three or five solutions to eq. (24). In 
particular, there is always a decaying solution, which we indicate with L. Its asymptotic behaviour is h{p) ^ (py/p) . 
This solution corresponds to a spacetime which is asymptotically fiat. Additionally, there are two or four solutions 
to eq. (24) which tend to a finite, nonzero value as p — ?> +oo. We name these solutions with C+, C_, Pi and P2 
(details about this denomination are given in [17]). Their asymptotic behaviour is h(p) — C, with C a constant. 
These solutions correspond to spacetimes which are asymptotically non-fiat. Interestingly, the leading term in the 
gravitational potentials scales as p^ for large radii, the same scaling which we find in a de Sitter spacetime. It is 
worthwhile to point out that, since we are working on scales below the Compton wavelength of the gravitational 
field, asymptotically non-flat does not really mean the real behaviour at infinity. To understand the true asymptotic 
behaviour of this solution, one should solve the complete, non-truncated equations. 

In a neighbourhood of p — >■ 0+ there are either one or three solutions to eq. (24) . For /? > there are exactly three 
inner solutions, while for /? < there is only one inner solution. In particular, there is always a diverging solution, 
which we denote by D. Its leading behaviour is h{p) ^ — \/2/l3 {py/p). This solution exists for both /3 > and 
/? < 0, with opposite signs for each case. Using this solution in eqs. (22)-(23), one realises that the h^ term cancels 
the GM/p term, so the gravitational field is self-shielded and does not diverge as p — ^ 0+. This solution is in strong 
disagreement with gravitational observations. For /3 > 0, there are two additional solutions to eq. (24), which tend 
to a finite, non-zero value as p ^> 0"^. We indicate these solutions by F+ and F_ , and their leading behaviour is 
h{p) = ±(3/3)~^/^. Notice that for /3 < there are no solutions to eq. (24) which tend to a finite value as p — ?> 0+. 
The expressions (22)- (23) for the gravitational potentials imply that the metric associated to these solutions (F+ and 
F_) approximate the linearised Schwarzschild metric as p — > 0^. 

From the behaviour of the inner solutions, one concludes that only in the /3 > part of the phase space solutions 
may exhibit the Vainshtein mechanism [18], but not necessarily for all values of a [17]. The phase space diagram 
which displays our results about solution matching is given in figure 2. We discuss separately the /3 > and /3 < 
part of the phase space, and refer to the figure for the numbering of the regions. The notation I o A means that 
there is matching between the inner solution I and the asymptotic solution A. 





FIG. 2: Phase space diagram in (a, /?) for the solutions to the quintic equation (24) in h, where the different regions show 
different matching of inner solutions to asymptotic ones. The lines splitting the regions are half parabolas (/3 oc a^, with a > 
or a < 0) due to rescaling symmetry of eq. (24) . 



/?<0 

In this part of the phase space, there is only one inner solution, D, so there can be at most one global solution to 

(24). There are three distinct regions which differ in the way the matching works (see [17] for details): 

- region 1: D o C+. The boundaries of this region are the line {3 — ior a < and the parabola /3 — C12 a^ for 



a > 0, where C12 is the negative root of the equation —4 — 8y 



1076y^-K2883y'' = 0. 



- region 2: No matching. The boundaries of this region are the parabola /3 = ci2 a^ and the parabola /3 = c_ a^, 
where c_ is the only real root of the equation 8 + 48 y — 435 y^ + 676 y^ — 0. 

- region 3: D f^ P2. 

/3>0 

In this part of the phase space, there are three inner solutions, D, F+ and F_, so there can be at most three global 

solutions to eq. (24). There are six distinct regions with different matching properties: 

- region 4: F_ f^ L , D f^ C_. The boundaries of this region are the parabola (3 — C45 a^, where C45 = 1/12, 
and the line (3 = 0. 

- region 5: F_|_ ^^ C_|_ , F_ <-^ L, D o C_. The boundaries of this region are the parabola (3 = C45 a^ for a > 
and the parabola f3 = c^a o? for a < 0, where C56 = (5 + -^13)724. 

- region 6: D o C_ , F+ o C+. The boundaries of this region are the parabolas /3 = csg a^ and (3 = cgya^, 
where cgy is the positive root of the equation —4 — 8 j/ + 88 y^ — 1076 y^ + 2883 y"^ = 0. 

- region 7: F+ ^ C+. The boundaries of this region are the parabola /3 = cg? a^ and the parabola (3 = 0+0^, 
where c+ = 1/4. 

- region 8: F+ o C+ , F_ o Pi , D f^- P2. The boundaries of this region are the parabolas /? = c+ a^ and 
13 = csg c?, where cgg = (5 - \/l3)/24. 

- region 9: F_ f^ Pi , D f^ P2. 

We note that the decaying solution L never connects to the diverging configuration D, so we can not have a spacetime 
which is asymptotically flat and exhibits the self-shielding of the gravitational field at the origin. On the other hand, 
finite non-zero asymptotic solutions (C± or Pi, 2) can connect to both finite and diverging inner solutions. Therefore, 
one can have an asymptotically non-flat spacetime which presents self-shielding at the origin, or an asymptotically 
non-flat spacetime which tends to Schwarzschild spacetime for small radii. More precisely, for /? < there are only 
solutions displaying the self-shielding of the gravitational field, apart from region 2 where there are no global solutions. 
Therefore the Vainshtein mechanism never works for /3 < 0. In contrast, for /3 > all three kinds of global solutions 
are present. Solutions with asymptotic flatness and the Vainshtein mechanism are present in regions 4 and 5, while 
solutions which are asymptotically non-flat and exhibit the Vainshtein mechanism do exist in all (/3 > 0) regions but 
region 4. Finally, solutions which display the self-shielding of the gravitational field are present in all (J3 > 0) regions 
but region 7. 

For the sake of clearness, we show one representative plot with the numerical matching solution between the inner 
and the asymptotic solutions. We consider solutions which recover the Schwarzschild solution near the origin, and 
which are asymptotically flat (F_ o L in Figure 3-left), or non-flat (F_|_ ^^ C+ in Figure 3-right). Finally, it is 
essential to decide whether these vacuum solutions are indeed consistent with matter sources, as it was done for the 
/3 = case in [16]. 

The Vainshtein mechanism has been studied intensively in the context of the Dvali-Gabadadze-Porrati braneworld 
model [24] and Galileon models [19]. Especially, it was shown that the most general second order scalar tensor theory 
described by the Hordenski action [20] leads to the same field equations as massive gravity [21]. 

B. Branch II: exact solutions 

As we learned in the previous section, an essential property of this theory of massive gravity is the strong coupling 
phenomenon occurring in the proximity of a source. On the other hand, the graviton mass induces non-linearities in 
the behavior of long wave length gravitons, responsible for the emergence of the second branch of solutions that we 
are going to study in this section. In an appropriate gauge, these solutions are asymptotically de Sitter or Anti-de 
Sitter, depending on the choice of parameters. 

We start with the unitary gauge (tt^ = 0) and allow for arbitrary couplings as and a^, while from now on we set, 
for simplicity, the bare cosmological constant to vanish (see [12, 15] for a more complete discussion including a bare 
cosmological constant). We choose the static Ansatz of eq. (12) for the metric, and we focus on the second branch 
of solutions for the constraint equation (13): c{r) = CQr. Then, the exact solution of field equations is given by 



nVn',„ 



f/fa. 



P/P„ 




FIG. 3: Numerical solutions for F- o L (left) and F+ ^ C+ (right), with (q, /3) = (0 , 0.1). foR = naR = -2GM/p. Solutions 
represent the GR gravitational potentials, and clearly present the vDVZ discontinuity with its resolution via the Vainshtein 
mechanism. Solutions asymptote a flat spacetime but with a different 7 — f/n than the Schwarzschild solution in GR (left) 
and a non-flat spacetime (right). 



[12, 15, 22] 



c{r) = Cor, 

b{rf = bo + 

a{rf+b{rf = Qo , 

d^ir)+airfbirf = Aq . 



h 
r 



62 r^ 



(30) 

Moreover, the equations of motion fix the constant parameters 60 1 ^2 , co , Qo , leaving the values of 61 , Aq free (although 
their sizes must be contained within certain intervals). Notice that in General Relativity, diffeomorphism invariance 
allows one to choose the function c(r) to be c{r) = r, so that Cq — 1. In this theory of massive gravity, after having 
fixed the gauge, this choice is no longer possible and the equations of motion determine cq. One finds 



Co 



1 + 60:3 + 12^4 ± yjl + 3a3 + 9a§ - 12^4 
3(l + 3a3 + 4a4) 



(31) 



for non a3 ^ — 4a4, and 



2 /I 



^"=3 



1 



12a4\ 
- 8q;4 / 



for as = — 4a4, which in particular includes the case a^ = a4=0. After plugging the metric components (30) in the 
remaining Einstein equations, one can find the values for the other parameters. The corresponding general expressions 
are quite lengthy, and for this reason we relegate them to Appendix A. As a concrete, simple example, in the main 
text we work out the special case a^ = —4 014, where the parameters are 



Co = 



1 



12a4\ 



1 — 8a4 
m^Ao 



4(12^4- 1) 



&o = 
Qo 



Ao 



16(1 - 12a4)^ + 81(1 - 8a4)^Ao 
36[l+4a4(-5 + 24a4)]^ 



(32) 



The previous solution is valid for 04 in the ranges 04 < 1/12 and 04 > 1/8. We find that 61 and Aq are arbitrary; this 
vacuum solution is then characterized by two integration constants. The resulting metric coefficients can be rewritten 
in the following, easier-to-handle form: 

1 - 8^4 y 



a{rf 

b{rf 

d{r) 



9 A 



1- 12a4/ 



9 A 

9Ao 



1 



8a4 Y 



12a4/ 



80:4 



1 - 12a4 



[p{r) +7 + 1] 



[l-p(r)] 



\/p{r) {p{r) + 7) 



c{r) = - 



1-I2a4\ 

r 

1 - 8^4 / 



(33) 



with (fi ~ —bi/bo) 



^^^ - ~ 9(1- 8^4)2 ' ^ - 81Ao ^ 1 - 8^4 y ^ ' 



r 



In order to have a consistent solution, we must demand that the argument of the square root appearing in the 
expression for d{r), Eq. (33), is positive. A sufficient condition to ensure this is that fJ.> 0, and 

r— 4/1- 12a4\^ , , 

The metric might be rewritten in a more transparent diagonal form, by means of a coordinate transformation. In 
absence of diffeomorphism invariance, any coordinate transformation of time t forces us to leave the unitary gauge, 
and to switch on a non-trivial profile for the Stuckelberg field tt^ of the form n^ = (7ro(r), 0,0, 0). One finds that 
then the metric can be rewritten in a diagonal form, as 

ds^ = -b{rfdt^ + d{rfdr^ + c{rf dil^ , (36) 

while the equations of motion for the fields involved are solved by 

,M^.l(i^)^ 1 , ,;m.-V^SM^, ,37, 

9 \ 1 — 8a4 / 1 — p(r) I — p(r) 

with b{r), c{r) and p{r) being the same as in Eq. (33). If one then makes a further time rescaling 

^ ^ 4(l-12a4)^ ^ ^33^ 



9Ay^ (1 - 8^4)^ 

the resulting metric acquires a manifestly de Sitter-Schwarzschild, or Anti-de Sitter-Schwarzschild form. The choice 
between these two possibilities depends on whether a^ is smaller or larger than 1/12, as can be seen inspecting the 
function p{r) in Eq. (34). On the other hand, we should point out that this time- rescaling cannot be performed, 
without further introducing a time dependent contribution to ttq. As expected, the metric in Eq. (36) can also be 
obtained by making the following transformation of the time coordinate dt = dt + -K^dr to the original metric (12). 
This produces a non-zero time component for tt^, that does not vanish even in the limit m — > 0. 

To summarize so far, we found vacuum solutions in this theory that are asymptotically de Sitter or Anti-de Sitter, 
depending on the choice of the parameters. 

Figure 4 shows the allowed parameters a^ and a^ for the existence of these asymptotically dS or AdS solutions. 
Further solutions and studies on black holes in this massive gravity theory can be found in [23]. 

IV. COSMOLOGICAL ACCELERATION 
A. Self-accelerating solution 

One of the interesting features of massive gravity is self- acceleration. The self- accelerating solution was originally 
found in the DGP braneworld model where the acceleration of the Universe can be realised without introducing the 
cosmological constant [25]. However the self-accelerating solution in the DGP model suffers from a ghost instability 
[26]. 

The first complete self-accelerating solution in the A3 massive gravity theory (6) was reported in [12, 15] (the self- 
accelerating solution in the decoupling limit was first obtained in [27]). This configuration describes an accelerating 
cosmological universe in the vacuum, in which the rate of acceleration is controlled by the size of the graviton mass. 
The solution is a coordinate transformation of the exact solution (36)- (37), after having performed the time rescaling 
(38) (one can use (A1)-(A3) for more general values of as and 04). Let us review how to construct it for the simple 
case of 0:3 = — 4q!4. In the case of the asymptotically de Sitter solution (36)-(37) with 04 < 1/12 and /i = 0, the 
metric can also be written in a time dependent form, at the price of switching on additional components of tt^ . After 
dubbing 



m = 



m 



2 



(1 - 12^4) 
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FIG. 4: Allowed region of parameter space where solutions, which asymptote dS (vertical lines in light or red color) or AdS 
(horizontal lines in dark or blue color), exist. Left (right) plot is for the negative (positive) branch of Eq. (31). The brane 
solution is given by the special choice of for as = —1/3 and Q4 = 1/12 [12], and the solid lines correspond to the line Q3 = — 4q4. 
We have set A = 0. 



we can make the following coordinate transformation t ~ Ft (r, p) and r ~ Ft (r, p) with 

/sinh(^) 



Ft{T,p) 



-4^ (i^l^) arctanh 



3 / 1 — 8a4 \ 



a4 

mr 12 



Vcosh(^)- 



m^ p^ 



t/2 



Then, the metric becomes that of flat slicing of de Sitter 



ds' 



-dT^ + e"" 



{dp^ + p^dfl^) 



where the Hubble parameter is given by 



The Stiickelberg fields tt^ are now given by 



m 



H ^ - ^ - 

2 2(1 -1204)5 



(7ro + r-Ft,p- Fr,0,0)^ 



(39) 

(40) 

(41) 
(42) 

(43) 



Interestingly, the value of the Hubble parameter is ruled by the mass of the graviton: we have a self-accelerating 
solution, in which the smallness of the acceleration rate can be associated with the smallness of the graviton mass. 

So far, we used the general static solutions of the previous section to determine time-dependent self-accelerating 
configurations via suitable coordinate transformations. However, one can also follow another approach, and try to 
directly find time-dependent, self-accelerating configurations in the theory of A3 massive gravity without relying on 
the unitary gauge. The hope, following this second route, is to determine additional cosmological configurations for 
this theory. Starting from the Ansatz we wrote in eq. (10), more general cosmological solutions can be obtained by 
focussing on the Ansatz d{r,t) = and c{r,t) ~ ra{r,t) in (10) so that the metric becomes 



ds' 



~b{t, rfdt^ + a{t, rfidr^ + r^dfl^). 



(44) 



Wyman et al [28] showed that self-accelerating configurations are characterised by the following profile for the function 
g{r,t) characterizing the Stiickelberg field 0' (see eq (11)) 



7(i,r) ^CQ^a{t,r)r, 



(45) 
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where Cq is given by (31). The equation of motion for g evaluated on the solution (45) provides a constraint on the 
function / characterizing the Struckclberg field (jp 



where the P„ functions 



%/XPi' = (2coP2 - ^2 ) VK - ^0 , (46) 



Po{x) = -12-2a;(x-6) -12(x-l)(x-2)a3-24(a;- 1)^0:4, 

Pi(x) = 2(3 -2x)+6(a;-l)(x-3)a3 + 24(x- 1)^0:4, 

^2(2;) = -2 + 12(a;- 1)^3-24(2;- l)2a4. (47) 



are evaluated at x = Cq ^, P^ix) = dPn/dx. Moreover, 



X 



(i+^l) -(f+'-T)' »'^5;<^-''--'^''' <«' 



and /i = sgii{fg' — gf) where primes denote derivatives with respect to r and ovcrdots with respect to t. Using these 
equations, it is possible to show that the Einstein equations are given by 



2 _ 1 + 3a3 ± 2a5 2 



3(1 + 3a3 ± a^y 

where 

ag = 1 + 3q;3 + 9^3 - 12a4 . (50) 

Note that there are two branches of solutions. This approach leads to self- accelerating solutions where the Hubble 
parameter is determined by the mass of the graviton. 

There are many possible solutions for the function f{t,r) for given solutions of metric satisfying the Einstein 
equations (49). For example, in the simple case of 0:3 = a4 = 0, the configuration (41) and (43) is given by 

5=1, a = e-/^ fit, r) = -\e-^l\ + 1 (arctanh (^e^/^rn.) + arctanh ( _^ ^ ^^^^7^1^) ) ) ' ^^^^ 

Notice that this self-accelerating solution has a flat FRWL metric, however, one can also write it as an open or 
closed FRWL spacetime, with the price of changing the Stiickelberg fields accordingly. In all the FRWL frames, the 
Stiickelberg fields are inhomogeneous. In fact it was suggested that there was no FRW solution that keeps the FRW 
symmetry for the fiducial metric S(^ [29]. However, Gumrukcuoglu ct al found a special self-accelerating solution 
which represents an open universe where the fiducial metric S^,y respects the FRW symmetries of the physical metric 
9yi.v [30] . Their solution is given by 

^ = 1' ^- 1 T%na ^ /ft,^)^laoft) | + (^^)'/||; , ao(i) = sinh(mt/2). (52) 

1 — (mr)^/16 m 1 — (mr)^/16 

For this solution the fiducial metric preserves the FRW symmetry, E(^ = (9/4)diag((2ao/m)^, 1, 1, 1). The behaviour 
of perturbations around this particular self-accelerating solution is very different from the other solutions that break 
the FRW symmetry for the fiducial metric YJ^. At the linear order, scalar and vector perturbations have no kinetic 
terms, hence they are strongly coupled [31], which leads to non-linear instabilities [32, 33]. The absence of the scalar 
kinetic term originates from the special choice of the solution for / [34] that retains the FRW symmetry for the fiducial 
metric. This leads to an enhanced symmetry that eliminates the scalar perturbations [35]. In the rest of this review, 
we do not consider this class of self-accelerating solution and consider the case where the FRW symmetry is broken 
for the fiducial metric E^. However, we emphasise that the physical metric still retains the FRW symmetry in these 
solutions. 

Wyman et al [28] also showed that the ordinary Friedmann equation is obtained even if we add ordinary matter 
energy density Pm{t)- The matter only sees the effect of the mass terms as a cosmological constant with no direct 
coupling to the scalar fields on the exact solution. Cosmological solutions in massive gravity and its extension including 
de a Sitter fiducial metric and bigravity can be found in [36]. 
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B. Decoupling limit solutions and their instability 

Once we determined self- accelerating, de Sitter solutions in this model, it is crucial to study their stability: this is 
the subject of this section. In order to make the analysis manageable, we focus on a convenient limit of Lagrangian 
(6) which captures most of the dynamics of the helicity-0 and helicity-1 mode, but keeps the linear behaviour of the 
helicity-2 (tensor) mode [3]. The limit, called the decoupling limit, is defined as 

m ^ , Mpi ^ oo , A3 EEE rn^Mpi = fixed , (53) 

In order to obtain canonically normalized kinetic terms for the helicity 2 and helicity 1 modes, together with the 
relevant couplings for the helicity modes, when this limit is taken one needs to canonically normalise the fields in 
the following way 

h^,^ -^ Mpi hf,^ , Af, -^ mMpi A^, TT ^ m^Mpi tt, (54) 

where we have split the Stiickclberg fields tt^ into a scalar component tt and a divergenceless vector A^ in the usual 
way, namely 

7r^ = r;^'^(a,7r + A.). (55) 

In order to take the decoupling limit (53) of the self-accelerating solutions defined by (45), the solution has to be 
in the conformally flat frame, which is defined as follows 

a{t,r) = b{t,r) = ^^ = (1 + H^ir"^ - t')/4)"^ (56) 

r 

In this frame, all the known self-accelerating solutions lead to the same decoupling limit solution for the Stiickclberg 
fields, namely 



^0 ^ / l-f3a3 + a5 ^V 1 / (1 + 3^3 + 2^5) [(1 + 3q3 + a^Y - (2 + 3^3 + a^Y^^] ^^^^2 ^,2. 

V(2 + 3a3 + a5)yA;; ) 2^/ 3Ao(l + 3a3 + a5)*(2 + Sag + 05)2 '"'' ^"^ ' 

TT-- = - ^^.\ r + 0{m^), (57) 

1 + 3^3 + as 

with Qfs given by eq. (50). If we split tt^ into a scalar and vector piece as in (55), canonically normalise the fields as 
in (54), and take the decoupling limit (53), one gets [37, 38] 



2\ co^ (2 + 3^3 + "5)2(1 + 3^3 + 2^5) y 2co ' 

Ao = -^A3r\ (58) 

where cq is given by (31) (or (32) for as — —Aa^), and the Hubble parameter H by 

2 _ (1 + 3^3 ±2^5)^2 

3(l + 3«3±«5)2 ' ^ ' 

which is the generalisation of (42) for arbitrary a^ and a^. Moreover, there is a relation between Aq and Qo 

^ ^ (l + 3a3 + a5)^(l + 3a3 + 2a5) .g.^ 

° 3g2(2 + 3a3 + a5)2(l + 3a3 + a5)^ + (2 + 3a3 + a5)4(l + 3a3-h2a5)' ^ ' 

and in the case of AdS, there is an extra bound given by 

^2 ^ (2 + 3^3 + ^5)^(1 + 3^3 + 2^5) .„,x 

3(1 -I- 3a3 + 05) 

If one take the vector charge to zero Qo = (or equivalently if Aq = (1 + 3q;3 + a5)'*/(2 + 3q!3 + 05)'*), these solutions 
can be written in a simpler covariant way 

V = -^A3i72(x^a;^)r,^, 7: = ^^-Ia^x^^x^, A^" = 0. (62) 

z Z Co 
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Therefore, corrections of order rri^ in (57), which do not show in the decouphng hmit, are the main differences among 
solutions in the fuh theory. These solutions in the decoupling theory were also found in [27]. 

In this decoupling limit, the structure of the Lagrangian becomes much simpler, and the various self-accelerating 
configurations become the same. For these reasons, it is particularly convenient to study the dynamics of perturbations 
in this limit. If problems or instabilities arise in this limit, then are they unavoidably present also in the full theory 
outside the decoupling regime. Interestingly, it has been shown in [27] that, in the decoupling limit, the coupling 
between the scalar mode tt and the trace T of the energy momentum tensor vanishes around these self-accelerating 
configurations: hence, the coupling to matter is the same as in GR with no need to implement a Vainshtein mechanism. 
However, we have shown in [37, 38] that all these backgrounds present instabilities in the vector sector. We present 
here the main results concerning these instabilities, focussing on the case of a^ = a^ = 0. A generalisation to arbitrary 
values is straighforward and can be found in [37, 38] . We start by considering perturbations of the fields ft.^^ , A^ and 
TT which only depend on time and radial component. Namely 



^M-- = f^li^ + ^, 



fii'j 



^^ = A^ + i^, 



TT = TTq + TT, 



(63) 



where the background quantities (those with an index 0) are given by the self-accelerating solution (62). The La- 
grangian for the tensor and scalar perturbations (without further truncations) reads 



where (in units in which A3 = 1, that we adopt from now on) H^ — i and A^^ is given by 



(64) 



Xau 



u^^Ux. - nn^, + -r?^, [u' - n;;n;; 



(65) 



We can use the following field redefinition to decouple the helicity 2 from the helicity field: 



^^7/ - dfj_TTdi,TT . 



(66) 



Then the kinetic terms for tensor and scalar are diagonalized resulting, up to total derivatives, in 



(67) 



Let us emphasize that the previous Lagrangian contains terms which are quadratic on /i^^, but higher orders in the 
scalar field tt. The scalar field terms are the so called Galileon combinations. On the contrary, as mentioned before, 
the vector piece has an infinite number of interactions [37, 38]. For our purposes, it is enough to stop at fourth order 
in the fields, resulting in the Lagrangian 



C, 






(68) 



where F,, 



d^A, 



d^Afj^. As mentioned above, the kinetic terms for the vectors vanish; however, vectors become 
dynamical by coupling them with the scalar at third or higher order in fluctuations (this was already pointed out in 
[27, 37, 39]). Nevertheless, one should worry about higher derivatives in the equations of motion, since the previous 
Lagrangian contains contributions with two time derivatives in the scalar field tt. For systems coupling scalars with 
vectors, it is possible to find the combination that ensures that the equations of motion do not contain at all terms 
containing more than two time derivatives. It is a generalization of Galileon combinations which was explored in [40] 
and dubbed p-form Galileons. Up to fourth order in perturbations, the correct combination (without including higher 
derivatives in A^) is 



^p-form ^ ^^ [trntrF^ - 2trnF2] +60 jtrF^ trH^ - (trH) 



4trntrnF2 + AtTH^F'^ + 2trnFnF 



} (69) 



where uq and bo are arbitrary coefficients, and we have used the tr notation to simplify the index structure. The 
above third order action with the aforementioned properties was presented in [40] , while the fourth order one is as far 
as we know new. Comparing (68) with (69) we notice that while the third order action has the correct structure to 
avoid higher order time derivatives in the equation of motion, the fourth order Lagrangian docs not seem to satisfy 
this requirement. However, a suitable field redefinition allows to remove the fourth order term from the third order 
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contribution, leaving a healthy Lagrangian without higher derivative equations of motion (in agreement with the 
ghost-free statement of the theory) . 



On the other hand, although our scalar- vector Lagrangian (68) does not lead to a propagation of a sixth ghost mode, 
it does generally lead to a ghost-like instability around self-accelerating configurations, in which the ghost is one of 
the available vector modes. It has be shown in [37] that, when turning on a non-trivial profile for the background 
vector field, the corresponding Lagrangian for perturbations around the resulting configuration acquires kinetic terms 
for the vector with the wrong sign. Here, following [38], we instead directly point out the instability by analyzing 
the Hamiltonian associated with the Lagrangian obtained by combining the third order Lagrangian contained in (68) 
with the scalar kinetic term: 



£ 



third 



1 



3H^TTDTT--[n7rF^,F^"' 



6 



2d^,TTF>'PF''p] 



(70) 



where we have removed the hats over the field to simplify the notation. We choose, for simplicity, the gauge Aq = 0, 
diAi ~ 0, and by doing a standard (3 + l)-dccomposition, the previous Lagrangian reads 



C = -3H^ 



2tt Ai/\Ai + Att A^ - TT a A^A 



tj ^i^j 



(71) 



where the dots represent the terms without time derivatives, which we do not include since they do not play a role in 
the present discussion. The conjugate momenta to n and Ai are 



n^ = -6H^ 



1 



%m 



AAA 






7T A Aj + Att Ai - tt^ ^ Aj 



(72) 



In order to analyze the associated Hamiltonian, it is convenient to introduce the matrix Kij = Air Sij^T:,ij. li Hij = 0, 
then, we can easily invert the relations that define the conjugate momenta, and obtain the following Hamiltonian 



n 



iTi 



1 



12 7J2 l2i/2 

3 A.Ha, 
2 "" A, A A- 



H, 



27 H^ 



dH'^AiRA, 

A,. A A 



A^Ua^ 

A AA 



(73) 
(74) 



where the dots represent terms without momentum variables. The previous Hamiltonian is linear in H^; hence it 
is unbounded from below. Notice that this argument holds even in the limit in which H^ vanishes. In conclusion, 
perturbations of the background self-accelerating solution, along the direction of scalar fluctuations such that k^ — 0, 
admit unstable directions along which the system falls towards regions where the energy is unbounded from below. 
Similar conclusions hold for more generic Kij. Let us, for example, consider a Kij that is non- vanishing, and invertible. 
Then, after straightforward manipulations, one can show that the Hamiltonian can be written as 



■ n 



1 



gm 



AAk/AA, 



(n. 



AA.k^/Ua^ 



n^.^^'n^. 



(75) 



where, again, the dots represent terms without momentum variables. It is not difficult to see that there are many 
unstable directions associated with this Hamiltonian. For example, make a choice for the vector AAi so that the 
scalar combination C = AAiK~- AAj is non-vanishing and has a given sign. For definiteness, the magnitude of AAi 
is chosen such that the denominator of the first term has the same sign of C. Accordingly, choose the vector H^i^ such 
that H^^K^-^n^^ has the same sign of C (for example, choose it in the same direction of the A^l^). Then, by choosing a 
suitable magnitude for H^r , it is possible to make one of the two terms in the previous Hamiltonian arbitrarily negative 
- hence the Hamiltonian is unbounded from below. Other cases, such as the case in which k~- is non-vanishing but 
not invertible, can be treated in a similar way. Furthermore, while here we focussed on the case a^ = a^ — it is 
straightforward to extend this analysis to the more general case, obtaining the same conclusion (see [38] for details). 
To summarize, one generically expects instabilities around the self-accelerating solutions discussed so far: there 
are many directions in the moduli space of fluctuations along which the energy is unbounded from below, and 
towards which the system can be driven into dangerous regions. On the other hand, to close this section with a 
positive perspective, it might very well be that suitable deformations of known solutions (or even completely new 
configurations) exist that, renouncing to the symmetries imposed on the Ansatze considered so far, do not present the 
problems discussed above. Very recently, a proposal in this direction has been pushed forward in [33], that consider 
the possibility of breaking the isotropy of three dimensional spatial slices to find stable configurations. Still much 
work is needed to clarify this subject and analyze phenomenological consequences of these solutions. 
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V. FUTURE DIRECTIONS 

Massive gravity is a good theoretical laboratory to study niodifieations of General Relativity with interesting 
phenomenological consequences. Non-linear self-interactions of massive gravity in proximity of a source manage to 
mimic the predictions of linearised General Relativity, hence agreeing with solar-system precision measurements. 
Moreover, massive gravity offers a concrete set-up for studying models of dark energy in modified gravity scenarios. 
Indeed, at large distances gravity is modified with respect to GR, and the theory admits cosmological accelerating 
solutions in the vacuum in which the size of acceleration depends on the graviton mass. Dark energy models built in 
this way have the opportunity to be technically natural in the 't Hooft sense: in the limit of graviton mass going to 
zero one gains a symmetry, by recovering the full diffeomorphism invariance of GR. Consequently, any corrections to 
the size of dark energy must be proportional to the (tiny) graviton mass itself. 

Hence, non-linear effects play a crucial role for characterizing phenomenological consequences of massive gravity. 
Motivated by this fact, the analysis of exact solutions of the equations of motion, obtained by imposing appropriate 
symmetries (spherical symmetry for static space-times, or homogeneity and isotropy for cosmological set-ups), make 
manifest, in idealized but representative situations, how the non-linear dynamics of the graviton degrees of freedom 
respond to the presence of a source, or, at very large scales, to the graviton mass itself. This has been the argument 
of this article, in which we reviewed our works on these topics. 

Much interesting work is left for the future: our results can be extended in various directions that will improve 
our understanding of massive gravity and, in general, of consistent infrared modifications of General Relativity. From 
one side, it would be interesting to find new stationary configurations renouncing to spherical symmetry, to test 
analytically the effectiveness of Vainshtein mechanism when spherical symmetry is broken. As a concrete example, it 
would be interesting to find analogues of the Kerr geometry in this scenario, in which frame dragging effects can be 
quantitatively analyzed. Also, it would be interesting to determine cosmological configurations that break the isotropy 
or homogeneity of the cosmological solutions analyzed until now. Indeed, working in a suitable decoupling limit, we 
have shown that the cosmological self-accelerating configurations studied so far are characterized by instabilities in 
the vector sector. Given the recent results of [33], we speculate that these instabilities can be possibly avoided 
by renouncing to some of the symmetries that characterize the solutions (for example isotropy of the three spatial 
directions). It would be interesting to determine stable self-accelerating backgrounds following this route, and study 
their consequences for what respect the dynamics of cosmological fluctuations. We hope to be able to develp all these 
questions in our future work. 
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Appendix A: General exact solution 

From the general Lagrangian (6), and using the non-diagonal ansatz (12) together with Einstein equations (7), 
one can show that there are two branches of solutions for non-vanishing a^ and 0:4 as it was done in Section III. 
Here we only consider the branch with a non-diagonal metric, where analytic solutions can be found. Since the 
combination y/l + Saa + 9a^ — \2a^ is always present in the solution of this branch (see (31)), it is convenient to 
map the (as, 04) parameters into (03, as), where a| = 1 + iaj, + 9q!§ — Yla^. In this new set of parameters, the 
combination d(r) Gtt + biry^ Gtr = 0, fixes c{r) as a function of r in the following way 

f \ (1+3^3 + ^5) ..,^ 

c(r) — Cor — — -r. (Al) 



(2 -f 3^3 + as) 



The rest of Einstein equations give 



Ao,, __^ ^ 2_^0, 



2 v-^ f;, "-v 7 2 



6(r)2 = -^(l-p), a(r)2 = -il(p + 7+l), d{r) ^ ^ ^^ - a{rfh(T)\ (A2) 

where 



/z (1 + 3a3 + 2a5) ,2 , _ (l + 3a3 + a5)^ 

^- r + 3(2 + 3a3 + a5)2"' ^ ' '^^^ Ao(2 + 3a3 + as)^ ^^'^> 
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Just like in the a^ — a4 — (a^ ~ 1) case, there are two integration constants, /i and Ag, but in order to have a 
positive argument for the square root in d{r), Aq has to run from Ag = to A™"^ = Cg. If we focus on the massless 
case fi — only, the solution describes the static patch of the de Sitter or Anti-de Sitter spacetime. 
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